Observation of a correlated free four-neutron system

A long-standing question in nuclear physics is whether chargeless nuclear systems can exist. To our knowledge, only neutron stars represent near-pure neutron systems, where neutrons are squeezed together by the gravitational force to very high densities. The experimental search for isolated multi-neutron systems has been an ongoing quest for several decades1, with a particular focus on the four-neutron system called the tetraneutron, resulting in only a few indications of its existence so far2–4, leaving the tetraneutron an elusive nuclear system for six decades. Here we report on the observation of a resonance-like structure near threshold in the four-neutron system that is consistent with a quasi-bound tetraneutron state existing for a very short time. The measured energy and width of this state provide a key benchmark for our understanding of the nuclear force. The use of an experimental approach based on a knockout reaction at large momentum transfer with a radioactive high-energy 8He beam was key.

A long-standing question in nuclear physics is whether chargeless nuclear systems can exist. To our knowledge, only neutron stars represent near-pure neutron systems, where neutrons are squeezed together by the gravitational force to very high densities. The experimental search for isolated multi-neutron systems has been an ongoing quest for several decades 1 , with a particular focus on the four-neutron system called the tetraneutron, resulting in only a few indications of its existence so far [2][3][4] , leaving the tetraneutron an elusive nuclear system for six decades. Here we report on the observation of a resonance-like structure near threshold in the four-neutron system that is consistent with a quasi-bound tetraneutron state existing for a very short time. The measured energy and width of this state provide a key benchmark for our understanding of the nuclear force. The use of an experimental approach based on a knockout reaction at large momentum transfer with a radioactive high-energy 8 He beam was key.
A neutron can be bound either in an atomic nucleus or in a neutron star. The free neutron has a lifetime of just under 15 min and decays into a proton, electron and antineutrino. The system made of two neutrons, the dineutron, is known to be unbound by only about 100 keV. Whether multi-neutron systems can exist as weakly bound states or very short-lived unbound resonant states has been a long-standing question 1 . The next simplest system of three neutrons is less likely to exist owing to the odd number of nucleons and therefore weaker binding; yet, a recent calculation has suggested its existence 5 . Following these considerations, the four-neutron system, the tetraneutron, is an appropriate candidate to address this question. An overview of previous experiments and theoretical approaches is given in ref. 1 .
Numerous attempts have been made to find a hint for the existence of the tetraneutron as a bound or resonant state. Among these attempts, experiments were performed searching for possible bound tetraneutrons produced in uranium fission reactions (see, for example, ref. 6 ). Other attempts, sensitive to both bound and resonant states, used pion-induced double-charge-exchange (DCX) reactions, mainly the He(π , π ) 4 − + reaction (see, for example, ref. 7 ), as well as transfer reactions such as He(d, Li) 8 6 (ref. 8 ). None of the experiments yielded a positive signal.
Most of the past experiments were performed with stable nuclei. Towards the twenty-first century, with the development of radioactive-ion beam facilities, it became possible to use extremely neutron-rich nuclei in which one can expect an enhanced formation of a tetraneutron system. The first indication for a possible bound tetraneutron was reported in 2002 2 from a break-up reaction of 14 Be into 10 Be + 4 n. The result stimulated several theoretical studies, all agreeing on the same conclusion: a bound tetraneutron state cannot be obtained theoretically without significantly changing our understanding of the nuclear forces [9][10][11] . However, the possibility of the four-neutron system existing as a resonant quasi-bound state with a very short lifetime on the order of a few 10 −22 s, before decaying, has remained an open and challenging question. It was later found that the result reported in ref. 2 is also consistent with such a resonant state with the limit on its energy ≲ E 2 MeV r (ref. 3 ). A decade later, in 2016, an indication of a tetraneutron resonance was reported 4 . A DCX reaction was used, but in contrast to previous attempts, this time the reaction was induced by a high-energy 8 He radioactive beam. 8 He is the most neutron-rich bound isotope, and the 8 He( 4 He, 8 Be) reaction channel was investigated. The advantage of using a radioactive beam is the freedom of selecting the reaction partner in a so-called recoil-less production (without momentum transfer) of the four-neutron system. The energy of the state was found to be E r = 0.8 ± 1.4 MeV, and an upper limit on its width was estimated as Γ ≤ 2.6 MeV. However, owing to the large experimental uncertainty, the possibility of a bound state could not be excluded by this experiment.
In this work, we used the quasi-elastic knockout of an α-particle ( 4 He nucleus) from a high-energy 8 He projectile induced by a proton target to populate a possible tetraneutron state. The inverse-kinematics knockout reaction He(p, p He) 8 4 at large momentum transfer is well suited because the 8 He nucleus has the pronounced cluster structure of an α-core ( 4 He) and four valence neutrons with small 4n centreof-mass motion, such that after the sudden removal of the α-particle, a rather localized four-neutron system with small relative energy between the neutrons is produced, which may have a large overlap with a tetraneutron state 12,13 . The chosen kinematics at large momentum transfer between the proton and the α-particle ensures that the four-neutron system will recoil only with the intrinsic momentum of the 4 He core in the 8 He rest frame, without any further momentum transfer, thus allowing the recoil-less production. Furthermore, final-state interactions between the four neutrons and the charged particles are also minimized owing to the large momentum transfer, separating charged reaction partners from the neutron spectators in momentum space (Fig. 1).
The experiment took place at the Radioactive Ion Beam Factory operated by the RIKEN Nishina Center and the Center for Nuclear Study, University of Tokyo, using the Superconducting Analyzer for Multi-particles from Radio Isotope Beams (SAMURAI) 14 . A primary beam of 18 O was directed onto a beryllium production target producing a cocktail of radioactive nuclei from fragmentation. The secondary 8 He beam was separated using the BigRIPS fragment separator and transported with an energy of 156 MeV per nucleon to a 5-cm-thick liquid-hydrogen target 15 located at the SAMURAI spectrometer (Fig. 2).
The incoming beam was measured upstream of the target on an event-by-event basis using scintillators for charge identification as well as momentum measurement, and two drift chambers for tracking (Extended Data Fig. 1).
The outgoing charged fragments (α-particle and proton) emerging from the quasi-elastic scattering were detected using a combination of detectors downstream of the target. Three planes of silicon-strip detectors, where each plane consists of two orthogonal layers enabling position measurements in both horizontal and vertical directions, served for tracking, energy-deposition measurement and reconstruction of the reaction vertex inside the target (Extended Data Figs. 2 and 3).
Behind the silicon planes, both charged fragments were bent through the magnetic field of the SAMURAI spectrometer, which was operated at a nominal magnetic field of 1.25 T in the centre of the magnet. The experiment was designed to detect an α-particle and a proton that emerge from quasi-elastic scattering close to 180° in the centre-of-mass frame ( Fig. 1). Under these kinematical conditions, their resulting outgoing momenta are very different from each other in the laboratory frame, as shown in Figs. 1 and 2. The knocked-out α-particle is slowed down from its initial momentum, that is, with the incoming beam velocity, to a momentum of about 330 MeV/c per nucleon after the reaction (where c is the speed of light). In contrast, the proton, which was at rest in the initial state, becomes the fastest particle in the reaction, gaining a typical momentum of about 860 MeV/c. At the focal plane, a drift chamber is used for tracking of the fragments after the magnet, and two scintillator walls located side by side, which cover a wide momentum range, are used for energy-deposition and time-of-flight measurements. The α-particle and proton are identified from a combination of their measured energy deposition, each in a different scintillator wall, and their position in the drift chamber (Extended Data Fig. 4). Their momenta are determined precisely from their reconstructed trajectories through the SAMURAI spectrometer.
As no additional momentum is transferred to the neutrons in the reaction, they continue moving with nearly beam velocity and can be detected, in principle, by the neutron detectors placed at a forward angle behind the SAMURAI spectrometer. The detection efficiency for neutrons is much lower than that for charged particles and decreases quickly as a function of the number of detected neutrons. The small p-4 He elastic cross-section at backwards centre-of-mass angles of less than 1 microbarn (ref. 16 ) resulted in the relatively low statistics of 422 events obtained for the He(p, p He) 8 4 reaction. These factors made it impossible to detect more than two neutrons in coincidence with the charged particles. Therefore, the neutron detection is not a part of the current study, aside from a consistency check (provided in Supplementary Information) of the near recoil-less production of the free neutrons.
The combined selection of event-by-event identification of incoming 8 He-beam particles in coincidence with the knocked-out α-particle and the scattered proton defines the He(p, p He) 8 4 channel. From a precise measurement of the momenta of the charged particles, the energy spectrum of the 4n system is reconstructed assuming energy and momentum conservation through the missing mass:  In this frame, the momentum of the proton balances that of the 4 He, P P = − p Article this notation, a bound 4 n system will appear at E 4n < 0 whereas a resonant state will appear at E 4n > 0. The missing momentum in equation (1)  , where the four-momenta P on the right-hand side of the equation are those of the incoming beam, target proton, knocked-out α-particle and scattered proton, respectively.
The He(p, p He) 6 4 knockout reaction was measured with almost exactly the same experimental conditions as for 8 He, except for some small differences in the energy of the incoming beam and the beam profile (Supplementary Table 2), and served as a benchmark for verifying the analysis and calibration procedures. In the case of 6 He, the 2n system is produced by the sudden removal of the 4 He core. The two-neutron relative-energy spectrum is expected to be well described by theory taking into account both the well established ground-state wavefunction and the final-state scattering wave of the two neutrons, predicting a low-energy peak around 100 keV. Similarly to the 8 He case, we define the missing mass (P P → He He 8 6 and m m 4 → 2 n n ). The measured missing-mass spectrum for 6 He is shown in the right panel of Fig. 3 together with the theoretical calculation 17 convoluted with the experimental acceptance and resolution (blue curve). The energy range shown represents the one covered by the experimental acceptance. The calculation is compared with the data by implementing it into an event generator for the quasi-elastic reaction, which uses the measured p-4 He differential elastic cross-section 16 as an input, as well as the measured internal momentum distribution of the α-particle in 6 He (ref. 18 ). The generated events are transported through the experimental set-up in Geant4 simulations to account for the experimental acceptance and detector resolutions. The excellent agreement of the simulated theoretical distribution with the measured spectrum confirms the analysis and the calibration for determining the missing mass. The missing-mass resolution obtained in the measurement is approximately 1 MeV sigma, and is almost constant over the measured energy range. The systematic uncertainty for the determination of the absolute energy was estimated from this measurement to be 0.4 MeV and that of the energy width to be 0.27 MeV (Methods). Also shown in the right panel of Fig. 3 (green curve) is a possible small background contribution coming from two-step process where 4 He is produced in a first step (see Methods and following discussion for 8 He). This background was estimated from the measured cross-section to contribute 1% of the total number of measured events.
The measured missing-mass spectrum of the four-neutron system from the He(p, p He) 8 4 reaction is shown in the left panel of Fig. 3.
Two components are observed: a well pronounced peak in the low-energy region with an energy around 2 MeV and a broad distribution at higher energies attributed to a non-resonant continuum response 13 , a direct four-body decay.
The shape of the non-resonant continuum spectrum of the four neutrons has been studied theoretically for the case of the four-neutron structure formed after the sudden removal of the α-core from 8 He (ref. 13 ). The creation of the system is investigated by introducing into the Schrödinger equation a source term that accounts for the reaction mechanism producing the four-body system, and that depends explicitly on the internal structure of the parent nucleus 8 He. The 8 He ground-state wavefunction (without final-sate interaction) was treated using the five-body ( He + 4n 4 ) cluster orbital shell model approximation (COSMA) 12 . The exact shape of the non-resonant continuum is sensitive to the hyperradius of the source, ρ sour an internal radius of the 4n system, described in the hyperspherical harmonics basis. A hyperradius of 5.6 fm is considered by the theory as the most realistic, as it reproduces the correct experimental radius of 8 He in the COSMA model. This results in a broad distribution centred around 30 MeV, in good agreement with the observed experimental spectrum.
We model the spectrum as follows: where a, b, and c are constants, f res is a Breit-Wigner function representing the possible resonance structure, and f con is the non-resonant continuum part presented above with the hyperradius as a parameter. The last term in equation (2), f bkg , represents possible background events coming from competing processes. Several processes were investigated and quantified (Methods), where the only non-negligible contribution found is from a two-step process involving 6 He ( 4 He) production: proton-induced break-up of 8 He into 6 He ( 4 He) followed by a p-4 He quasi-elastic scattering. The resulting energy distribution is broadened and shifted to lower energies compared with the pure 6 He case (right panel of Fig. 3) owing to the two-step process, which has been taken into account in the simulation of f bkg . This background was estimated from measured cross-sections to contribute 2.6% to the total number of measured events (Methods), which has been used to determine the normalization constant c.
The experimental spectrum was then fitted with the energydependent function given in equation (2), where the fit function was convoluted with the experimental response, taking into account (Extended Data Fig. 5). The result of the χ 2 minimization is presented by the solid blue curve in the left panel of Fig. 3, together with the individual contributions. The statistical significance of the peak structure is well beyond the 5σ level (Methods).
The probability of populating a four-neutron system in a resonant state after the sudden removal of the α-core in 8 He is determined by the overlap of the 4n wavefunction in the final state and the relative motion of the four neutrons in the 8 He initial state. This overlap integral defines the ratio of cross-sections for the population of the resonance and the non-resonant continuum. Unconvoluting with the acceptance of the set-up, following the energy dependence of equation (2), we extract a probability of P r = 18.7 ± 2.3%. For comparison, the relative motion of the four neutrons studied in the COSMA model 12,13 yields a probability of about 30%. This value is obtained by considering the hyperradius of 5.6 fm, whereas the resulting value from the fit to the experimental data is 5.0 ± 0.1 fm, which would yield a smaller probability to populate the resonant state.
Assuming a resonant state, its energy and width as determined from the fit are  4 (open red symbol). The energy of the resonance is in agreement within the uncertainty, despite the fact that different reactions were used to probe the 4n system, and is also in agreement with the upper limit given in ref. 3 .
From the theory side, there is no consensus among the different theories and their predictions are partly contradictory, although, there is a general agreement that a bound 4 n does not exist. In 2003, Pieper 10 studied this possibility using Green's function Monte Carlo calculations. His conclusion was that the existence of a bound 4 n state has to be excluded, unless nuclear forces are drastically modified. However, his calculation suggested that a possible resonance might exist near 2 MeV, but in such a case it must be very broad.
Using a similar approach, the quantum Monte Carlo (QMC) framework based on two-body and three-body chiral interactions was used to calculate the energy of a 4 n resonance 5 . The result supports the existence of a resonant state with an energy of 2.1(2) MeV, while no prediction has been made for its width (blue band). An extended no-core shell model (NCSM) approach using a harmonic-oscillator representation predicts different resonant states (full stars) including their corresponding widths 19 (see also ref. 20 cited in ref. 20 ). Calculations have been performed also in the framework of the no-core Gamow shell model (NCGSM) 21 . These resulted in E r ≈ 7 MeV and Γ ≈ 3.7 MeV (cross), where the conclusion was in fact that the energy of a 4 n resonance might be compatible with the experimental value of ref. 4 , albeit with a significantly  20 (cross) 21 , where the blue arrow indicates that the width is predicted to be larger than 3.7 MeV, and QMC calculations 5 (band). Whether this observation of a low-energy peak is attributed to a four-neutron resonant state or to other correlations between the neutrons in the final state, needs to be clarified by ab initio theories. Left: missing-mass spectrum of the four-neutron system extracted from the He(p, p He) 8 4 reaction. The different curves represent a Breit-Wigner resonance (pink), a non-resonant continuum (dashed blue), the background from two-step processes (green) and the total sum (solid blue). Right: missing-mass spectrum of the two-neutron system extracted from the He(p, p He) 6 4 reaction. The blue curve represents the theoretical calculation 17 convoluted with the experimental acceptance and resolution, and the green curve represents the background from the two-step reaction.
Article larger width. As pointed out in a later study 20 , these calculations were incomplete, as they were performed only in truncated model spaces or with unphysically overbinding interactions. The authors of this work 20 concluded that both the energy and the width of such a resonance are comparable with the experimental data (open star). At the same time, other calculations claim that to generate such a resonance, nuclear forces have to be significantly changed [22][23][24][25][26][27] , which would not be consistent with the present understanding. We note that some theories 26,27 predict a non-resonant low-energy enhancement of the density of states in the four-neutron spectrum. Whether such a prediction is consistent with our observed resonance-like feature cannot be currently ascertained, as the energy spectrum of the four-neutron system is not given. The drastically different predictions resulting from different theoretical approaches highlight the importance of the current firm experimental observation.
In conclusion, we have presented the experimental observation of a resonance-like structure consistent with a tetraneutron state near threshold after 60 years of experimental attempts to clarify the existence of this state. The use of a high-energy knockout reaction in inverse kinematics allowed a precise measurement. The use of a radioactive 8 He beam as the parent system and a direct, large momentum-transfer reaction opened up the opportunity to create the 4n system in a well defined one-step process and in a recoil-less undisturbed way. The optimized detection system enabled a precise determination of the final state and a high-resolution measurement. A well developed peak structure has been observed at an energy of 2.37 ± 0.38(stat.) ± 0.44(sys.) MeV with a striking statistical level. This is in agreement with the result of ref. 4 and the upper limit given in ref. 3 . Both the energy and the extracted width of Γ = 1.75 ± 0.22(stat.) ± 0.30(sys.) MeV are consistent with a tetraneutron state that is unbound with a corresponding lifetime of (3.8 ± 0.8) × 10 −22 s. Next-generation experiments using different reaction mechanisms and possibly detecting the four neutrons in coincidence will reveal more insights into the properties of the four-neutron system, including correlations among the neutrons. More elaborated ab initio nuclear theories accounting fully for the effect of the continuum are necessary to understand the observed low-energy peak and its origin.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-04827-6.

Methods
Energy E 2n of the 2n system The theoretical calculation used as an input for the E 2n distribution includes both the ground-state wavefunction and the final-state interaction between the two neutrons 17 . These were obtained from a three-body ( He + 2n 4 ) cluster model, which contains a phenomenological three-body force as well as the following local angular momentum l-dependent two-body interactions: a central s-wave interaction in the nn system, and central s-p-and d-wave interactions in the n-α system, as well as n-α spin-orbit interactions on the p-and d-waves. The ground-state wavefunction is computed in a three-body model using the computer code FaCE 28 , an acronym for Faddeev with Core Excitation, whereas the final-state interaction is taken into account by using the nn transition t-matrix approach. The calculated distribution is shown in Extended Data Fig. 6 (purple). For comparison, we also show the calculated distribution without the nn final-state interaction (green), which reflects the nn motion in the source, the 6 He ground state. The distributions are normalized to the same maximum, such that only their shape is compared. It can be seen that the effect of the nn final-state interaction is very large. The ground-state distribution peaks at a larger energy with a far extending tail towards higher energies.

Systematic uncertainty of the missing-mass determination
As the theoretical calculation used as an input for the energy of the two-neutron system in the quasi-elastic He(p, p He) 6 4 reaction is considered to be accurate, the comparison of it with the measured spectrum is used to estimate the systematic uncertainty of the missing-mass measurement. Different offsets covering the range of ±0.5 MeV are added to the generated theoretical distribution and for each offset the goodness of the agreement with the data is quantified by calculating the χ 2 value. The generated distributions are convoluted with the energy response of the set-up, to account for the experimental acceptance and the detector resolutions. The minimum χ 2 obtained is located very close to zero, that is, very close to the original distribution (Extended Data Fig. 7). The range of deviations in the energy scale still showing reasonable agreement with the data ( χ ± 1 min 2 range) is interpreted as the systematic uncertainty for the absolute missing-mass determination. This results in an uncertainty with an average value of ±0.4 MeV. To determine the systematic uncertainty for the missing-mass width, the original theoretical distribution is used (that is, offset of zero). To check the sensitivity to the resolution, different resolutions are applied in the range σ = 1± 0.5 MeV. The generated distributions are convoluted with the acceptance of the set-up and are smeared using a Gaussian function with a width σ. Similarly as described above, the χ 2 values are calculated and the χ ± 1 min 2 range is interpreted as the systematic uncertainty (Extended Data Fig. 7). This results in an uncertainty with an average value of ±0.27 MeV. We note that for both absolute value and resolution, the minimum χ 2 values fall very close to the expected values (which could be statistically by accident). As the given systematic uncertainties are determined solely by the statistics of the 6 He data, we consider the given uncertainties as conservative estimates.
An additional systematic uncertainty from the background function in equation (2), on both the energy and the width, is added quadratically. The fitting procedure is repeated three times, with the green curve presented in Fig. 3 and with the 1σ lower and upper limits (Supplementary Information). The values for the energy and the width are taken as the average between the different fits, and the systematic uncertainty is taken as the difference between them. This results in an uncertainty of 0.18 MeV and 0.14 MeV for the energy and width, respectively.
Energy E 4n of the 4n system The relative-energy spectrum of four neutrons produced by sudden removal of the α-particle from 8 He as a source was studied theoretically in ref. 13 . The final-state energy spectrum depends on the intrinsic 4n relative motion in 8 He, that is, the 8 He ground-state wavefunction as the source (without final-state interaction), and the final scattering state of the four isolated neutrons, that is, 4n final-state interaction. The wavefunction used in ref. 13 for 8 He is based on the COSMA model 12 , and the source term contains the Fourier transform of the overlap between 8 He and the α-particle wavefunctions. The method of hyperspherical harmonics 29 was used to solve the equations of the model. This method is based on the link between the hyperspherical-function method and the oscillator NCSM and uses a Slater determinant representation of the hyperspherical harmonics functions. The resulting source function is described by the internal variables of the four-neutron system, hyperradius, hyperangle, hyperangular momentum and the 4 He-4n relative motion. Calculations showed that only hyperangular momentum K = 2 significantly populates the 4n non-resonant continuum, reflecting the motion of the four neutrons in 8 He. The resulting energy spectrum in ref. 13 is presented only up to 20 MeV. It was also used in ref. 4 to estimate the non-resonant background and can be modelled analytically as where m N is the nucleon mass, ρ sour is the hyperradius, ħ is the reduced Planck constant and the factor 3.3 was estimated by matching to the calculated distribution in ref. 13 .

Background events from competing reactions
Events with a proton and an α-particle in the final state might result from other competing reactions. However, owing to the unique kinematics of quasi-elastic p-4 He scattering at large centre-of-mass angles, other direct reactions are excluded. Examples are one-neutron knockout, He(p, pn) He, 8 7 or 6 He knockout, He(p, p He*)2n. 8 6 In the former case, the momentum of the proton is too small to be selected by the experiment, where in the latter, it is too high at the angular range covered by the set-up, owing to the larger mass of 6 He relative to 4 He Therefore, the only possible contribution of background events to our measured p-4 He events can come from secondary reactions. Below we list the possible reactions and their expected contributions.
He(p, pn) He 8 7 . A first interaction is single-neutron knockout from 8 He leading to break-up into 6 He and two neutrons. This can be followed by a second interaction, quasi-elastic p-4 He scattering at backward centre-of-mass angles. These processes will result in events with a missing-mass distribution similar to that of the one-step He(p, p He) 6 4 reaction but broadened and shifted to smaller energies. This is due to the difference in the separation energies of two and four neutrons in 6 He versus 8 He of about 2 MeV, and due to the additional recoil (Supplementary Information). The single-neutron knockout cross-sections from 6,8 He were measured in an experiment performed within the same experimental campaign as the experiment presented here, at a slightly higher beam energy of 185 MeV per nucleon using a series of targets with different nuclear charges Z (ref. 30 ). For 6 He, the cross-section was extracted for the hydrogen target, whereas for 8 He, the lowest-Z target used was carbon. We therefore scale the measured cross-section for 6  . He(p, p) He* 8 8 . A first interaction is inelastic excitation of 8 He leading to a break-up into 6 He and two neutrons. These cross-sections were measured as well in ref. 30 , with scaling ratio of σ σ / =1.02 ± 0.13.
He,C He,C 8 6 Using σ = 11.0± 0.9 mb He,H 6 gives σ = 11.2± 1.7 mb inel . In total, the cross-section evaluated for the two processes described above equals σ = 71.1± 9.2 mb. tot Evaluating the relative contribution to the number of one-step p-4 He quasi-elastic reactions N reac at half of the target thickness (t) gives N N σ t / = (b) × (g cm )× bkg reac −2 0.6(Avogadro)/A = 0.071 × 0.175 × 0.6 = 0.75±0.10% , where A is the mass number of the target. These types of processes can also occur from 8 He break-up along beamline materials before the target, plastic scintillators, kapton windows and air. These contributions were estimated starting from the last plastic scintillators at the BigRIPS fragment separator up to the kapton window at the entrance to the target cell, using the measured cross-sections for hydrogen and carbon targets 30 . The produced 6 He has an angular spread in the transverse direction according to its intrinsic momentum, which has a width of σ c = 35 MeV/ He 6 (ref. 18 ). Assuming a central beam energy of 156 AMeV, its corresponding momentum is P c = 3.37 GeV/ He 6 , leading to an angular θ width of σ σ P = / = 0.035/3.37 = 10 mrad.
θ He He 6 6 The spread in the x and y directions evaluated at the entrance to the target equals σ d σ d = ⋅ = (mm) × 0.01(rad), where d is the distance from the material at which 6 He was produced to the target entrance. The target enclosure radius equals 20 mm. Therefore, depending on the distance of production point and target, some fraction of the produced 6 He is expected to not reach the target, which reduces its intensity. Using the measured distances to the target entrance, the estimated fraction of 6 He produced before the target and reaching the effective target volume is 0.36 ± 0.01%. Overall, the estimated contribution from these two processes is 1.11 ± 0.10%.
He(p, p He)2n 8 6 . A first interaction is 6 He knockout to its ground state. This process can contribute only when the 6 He is produced along the target, as, owing to its angular spread, it will not reach the target region in case of production in beamline materials. The cross-section for this process is not well known, and was studied previously only for a much higher energy 18 . However, it was shown to be compatible with the p-4 He elastic scattering cross-section. We therefore adapt the measured p-4 He total cross-section, which was measured at 156 MeV (ref. 16 ), σ = 91.8 mb, and use it as an upper limit. The relative contribution of these events is evaluated as 0.0918 × 0.175 × 0.6 = 0.96%.
He + p → He 8 4 . In the first interaction, 4 He is produced from fragmentation of 8 He. This can be, for example, a result of elastic p-4 He scattering or inelastic excitation of 8 He. The inclusive 4 He cross-section was measured at higher energy for He + C→ He 8 1 2 4 to be 95 ± 5 mb (ref. 31 ). We scale it by a factor of 1/2 to estimate the cross-section for the hydrogen target, leading to σ = 48 ± 5 mb. The relative contribution of these events is evaluated as 0.048 × 0.175 × 0.6 = 0.50± 0.05%. This process will result in events with a negative missing-mass distribution below −3.1 MeV, the binding energy of 4 He in 8 He that will extend to more negative values owing to the additional recoil (Supplementary Information).
He(p, p) He 8 8 . A first interaction is p- 8 He elastic scattering at backwards centre-of-mass angles. This will lead to a fast proton in the final state. In a second interaction, 4 He can be produced. However, in this case, the momentum of the proton is even larger than that resulting from direct 6 He knockout, owing to the larger mass; therefore, this reaction channel is excluded.
For the benchmark measurement with 6 He, only one reaction can contribute to a background, where 4 He is produced in the first interaction, He + p → He 6 4 . Similar to 8 He, the inclusive 4 He cross-section was measured at higher energy for He + C→ He 6 1 2 4 to be 189 ± 14 mb (ref. 31 ). We scale it by a factor of 1/2 to estimate the cross-section for the hydrogen target, leading to σ = 95 ± 7 mb. The relative contribution of these events is evaluated as 0.095 × 0.175 × 0.6 = 1.00± 0.07%. This process will result in events with a negative missing-mass distribution below −0.975 MeV, the binding energy of 4 He in 6 He that will extend to more negative values owing to the additional recoil (Supplementary Information).

Summary of background contributions
From the reactions listed above, we conclude that the only contributions from background processes come from reactions in which 6 He or 4 He were produced in a first step, and the overall contribution is evaluated as 2.6%. In the fitting of the missing-mass spectrum, f bgk is taken as the simulated distribution for the two-step processes described above. Each process is simulated individually, and adds weight to the total distribution according to its expected contribution. In addition, as the measured cross-sections adopted to estimate the background contributions were not measured directly for the reactions of interest, we consider the uncertainty as a variation by a factor of two on each one of the scaled cross-sections (Supplementary Information).
It is noted that the resulting background distribution can also explain the two events observed at the energy region of E −10 MeV < <−8 MeV. 4n These two events cannot result from a direct reaction, as the lower limit of a bound tetraneutron is −3.1 MeV, corresponding to the binding energy of 8 He against decay into He + 4n 4 . Therefore, these events can be attributed only to two-step reactions. Estimating the background function at −9 MeV, we expect 0.27 events in that energy region, leading to a 2σ significance of the measured two events, which is not considered as a statistically significant deviation from the expected background contribution in that energy region.
Finally, a major difference between one-step and two-step processes is expected in the reaction-vertex distribution of the proton and the α-particle (Extended Data Fig. 8). For a two-step process, an exponentially increasing yield of reaction products along the target is expected owing to the need to produce the 6 He ( 4 He) first. Such a distribution is not observed in the data, corroborating that the background contribution is indeed small. The estimated background contribution will be considered to evaluate the significance level of the observed peak.
To quantify the statistical significance of the resonance observed in Fig. 3, we evaluate the number of background events in the peak region defined as E Γ ± 2 r and E −2 MeV < < 6 MeV. 4n The number of measured events in the peak region amounts to N = 41 events. The number of background events equals the integral of the background and continuum functions (f bkg and f cont in equation (2)) in that region, resulting in an average number of eight events from the background and two events from the continuum, such that the total expected background amounts to ten events. Even though we allowed a variation of a factor of two for the background estimation, this does not change the result significantly. Hence, we conclude that the resonance structure is observed with high statistical significance, well above the 5σ observation threshold.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. Source data are provided with this paper.